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Abstract 

In  this  note  we  present  a  critical  review  of  the  some  of  the  positive 
features  as  well  as  some  of  the  shortcomings  of  the  generalized  sensitiv¬ 
ity  functions  (GSF)  of  Thomaseth-Cobelli  in  comparison  to  traditional 
sensitivity  functions  (TSF).  We  do  this  from  a  computational  perspec¬ 
tive  of  ordinary  least  squares  estimation  or  inverse  problems  using  two 
illustrative  examples:  the  Verhulst-Pearl  logistic  growth  model  and  a 
recently  developed  agricultural  production  network  model.  Because 
GSF  provide  information  on  the  relevance  of  data  measurements  for 
the  identification  of  certain  parameters  in  a  typical  parameter  estima¬ 
tion  problems,  we  argue  that  they  provide  the  basis  for  new  tools  for 
investigators  in  design  of  inverse  problem  studies. 
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1  Introduction 


The  interdisciplinary  character  of  scientific  research  today  leads  to  many 
instances  in  science,  engineering  and  finance  where  mathematical  models 
based  on  systems  where  either  algebraic  or  differential  equations  (or  both 
in  hybrid  models)  are  often  used.  In  recent  years,  because  of  an  increasing 
need  to  better  describe  social,  physical  or  biological  reality,  and  because  of  a 
tremendous  improvement  in  computing  technology,  the  models  employed  by 
scientists  have  became  more  complex  and  more  sophisticated,  wherein  the 
number  of  variables  and  parameters  being  utilized  have  increased  consider¬ 
ably.  This  naturally  has  motivated  new  questions  of  determining  the  relative 
importance  of  parameters,  the  effect  on  the  model  output  of  variation  in 
parameters,  the  uncertainty  in  the  model  results  due  to  the  uncertainty  of 
parameters,  etc.  Answers  are  not  obvious  for  large,  complex  models  and  the 
associated  problems  provide  significant  new  challenges  for  inverse  problem 
research  scientists.  While  the  traditional  literature  on  inverse  problems  has 
provided  a  wide  range  of  useful  theories  and  results  on  topics  from  well- 
posedness  and  regularity  techniques  to  computational  methodologies  (e.g., 
see  [23]  and  the  extensive  references  therein),  we  focus  here  on  concepts  re¬ 
lated  to  sensitivity  which  we  believe  can  aid  in  experimental  design  to  obtain 
data.  These  will  enhance  practical  aspects  of  parameter  estimation  as  well 
as  motivate  new  questions  in  basic  research.  We  do  this  in  the  context  of 
several  examples  in  which  we  introduce  and  illustrate  ideas. 

Sensitivity  analysis  is  an  ensemble  of  techniques  [24]  that  can  provide 
some  answers  to  questions  for  a  given  problem  of  interest,  yielding  a  much 
better  understanding  of  the  underlying  mathematical  model  with  a  resulting 
marked  improvement  in  the  results  obtained  using  the  models.  Tradition¬ 
ally,  sensitivity  analysis  referred  to  a  procedure  used  in  simulation  studies 
(direct  problems)  where  one  needed  to  evaluate  the  effects  of  parameter  vari¬ 
ations  on  the  time  course  of  model  outputs  and  to  identify  the  parameters 
or  the  initial  conditions  to  which  the  model  is  most/least  sensitive.  In  re¬ 
cent  years  however,  due  to  an  increasing  interest  in  incorporating  uncertainty 
into  models  and  in  ascertaining  the  sensitivity  of  parameter  estimates  with 
respect  to  data  measurements,  the  uses  of  sensitivity  have  broadened  sig¬ 
nificantly  [7,  26].  In  one  direction,  sensitivity  of  systems  with  probability 
measures  embedded  in  the  dynamics  (problems  involving  aggregate  dynam¬ 
ics)  have  become  important  in  applications  in  biology,  electromagnetics,  and 
hysteretic  and  polymeric  materials  (see  [7,  8]  and  the  references  therein).  On 
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the  other  hand,  investigators’  attention  has  also  recently  turned  to  the  sen¬ 
sitivity  of  the  solutions  to  inverse  problems  with  respect  to  data,  in  a  quest 
for  optimal  selection  of  data  measurements  in  experimental  design.  In  this 
paper  we  focus  on  this  latter  direction. 

As  part  of  model  validation  and  verification,  one  typically  needs  to  esti¬ 
mate  model  parameters  from  data  measurements,  and  a  related  question  of 
paramount  interest  is  related  to  sampling;  specifically,  at  which  time  points 
the  measurements  are  most  informative  in  the  estimation  of  a  given  param¬ 
eter.  Due  to  the  fact  that  in  practice  the  components  of  the  parameter  esti¬ 
mates  are  often  correlated,  traditional  sensitivity  functions  (TSF)  used  alone 
are  not  efficient  in  answering  this  question  because  TSF  do  not  take  into 
account  how  model  output  variations  affect  parameter  estimates  in  inverse 
problems.  In  an  effort  to  overcome  this  shortcoming,  Thomaseth  and  Cobclli 
[26]  recently  introduced  a  new  class  of  sensitivity  functions,  called  generalized 
sensitivity  functions  (GSF),  which  provide  information  on  the  relevance  of 
measurements  of  output  variables  of  a  system  for  the  identification  of  specific 
parameters.  For  a  given  set  of  time  observations,  Thomaseth  and  Cobelli  use 
theoretical  information  criteria  (the  Fisher  information  matrix)  to  establish 
a  relationship  between  the  monotonicity  of  the  GSF  curves  with  respect  to 
the  model  parameters  and  the  information  content  of  these  observations.  In 
this  paper  we  present  discussions  on  how  to  use  this  information  content 
tool  along  with  TSF  to  improve  the  parameter  estimates  in  inverse  problems 
and  therefore  to  further  validate  the  utility  of  these  new  functions  in  such 
problems.  It  is,  of  course,  intuitive  that  sampling  more  data  points  from  the 
region  indicated  by  the  GSF  to  be  the  “most  informative”  with  respect  to 
a  given  parameter  would  result  in  more  information  about  that  parameter, 
and  therefore  provide  more  accurate  estimates  for  it. 

We  present  a  critical  review  of  the  positive  features  as  well  as  shortcom¬ 
ings  of  the  GSF,  from  the  perspective  of  parameter  estimation  problems  via 
several  examples.  We  first  illustrate  ideas  with  the  classic  and  well-known 
logistic  growth  model  of  Verhulst-Pearl.  This  is  followed  by  discussion  in¬ 
volving  a  recently  developed  model  which  we  shall  refer  to  as  an  agricultural 
production  network  model.  Our  paper  is  organized  as  follows.  In  Section  2 
we  introduce  the  necessary  theoretical  framework  for  our  analysis  and  we  de¬ 
fine  the  traditional  (TSF)  and  the  generalized  (GSF)  sensitivity  functions.  In 
Sections  3  and  4  we  introduce  the  Verhulst-Pearl  logistic  growth  population 
model  and  the  agricultural  production  network  model,  respectively,  and  we 
discuss  numerical  simulations  carried  out  to  investigate  the  effectiveness  of 
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the  information  content  provided  by  the  GSF  to  improve  the  accuracy  of  the 
parameter  estimates  from  the  perspective  of  ordinary  least  square  problems 
with  noisy  data.  Finally,  in  Section  5  we  present  our  conclusions  and  remarks 
and  indicate  directions  for  future  work. 


2  Theoretical  Framework 

We  consider  a  parameter  estimation  problem  for  the  general  nonlinear  dy¬ 
namical  system 

xit)  =g(t,x(t),9)  ^ 

x(t0)  =  x0, 

with  discrete  time  observations 

Vj  —  f(tj,  0)  +  6j  —  Cx(tj,  9)  +  6j,  j  —  1, . . .  ,n  (2) 

where  x,  g  G  RN ,  /,  e3  G  MM  and  9  G  Mp.  The  matrix  C  is  an  M  x  N  matrix 
which  gives  the  observation  data  in  terms  of  the  components  of  the  state 
variable  x. 

We  assume  for  our  model  that  the  observation  errors  e3  are  independently 
identically  distributed  (i.i.d.),  with  zero  mean.  For  different  observation  co¬ 
ordinates  fi,  i  =  1, . . . ,  M,  we  have  different  variances  of  associated  with  the 
coordinates  of  the  errors  e3,  i.e. 

ej  ~  A/m(0,  V), 

where  V  =  diag(af, . . . ,  (72M).  We  also  make  the  standard  statistical  assump¬ 
tion  that  there  exists  a  “true”  value  parameter  90  such  that  the  data  set  {%•}, 
which  can  be  interpreted  as  a  realization  of  the  observation  process  Y  =  {Yj} 
at  the  discrete  time  points  tj,  j  =  1 , ,n,  has  the  form  in  equation  (2)  with 
9  =  90. 

We  use  an  ordinary  least  squares  (OLS)  approach  to  estimate  9q,  and  we 
seek  to  find  a  value  9n  that  minimizes  the  cost  functional 

n 

■w)  =  Y  (/(M)  -  yi)Ty~'  (/(v»)  -  vi)  ■  (3) 

3= 1 

Since  {yj}  is  a  realization  of  the  random  variable  set  {Yj},  the  estimate 
9n  we  obtain  by  minimizing  the  cost  functional  Jn  is  a  realization  of  some 
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random  variable  @n.  Therefore,  the  accuracy  of  our  parameter  estimates 
9n  ultimately  depends  on  the  statistical  properties  of  this  random  variable, 
and  in  order  to  qualitatively  analyze  the  estimates,  we  use  a  standard  error 
approach  [16,  19,  20,  25]. 

From  the  asymptotic  theory  of  statistical  analysis,  one  finds  that  as  n  — > 
oo,  0n  ~  Nv(9q,  S0)  is  a  good  approximation,  where  the  covariance  matrix 
S0  is  given  by 

E»=  .  (4) 

and  Dj{6 )  is  the  M  x  p  Jacobian  matrix  of  /  with  respect  to  9  at  tj ,  i.e., 


D‘(e)  ~  ~m~' 

These  are  quite  clearly  the  TSF  for  observations  or  measurement  outputs 
with  respect  to  the  parameters.  The  covariance  matrix  S0  is  used  in  formu¬ 
lating  the  standard  errors  for  our  estimates  9n\  these  are  given  by 

SEfc  =  V (So)fcfc,  k  =  1,2,  (5) 

Because  9q  in  (4)  is  unknown,  we  replace  it  by  9n  when  calculating  approx¬ 
imations  for  (5).  Moreover,  the  variances  af, ... ,  a2M  are  also  generally  un¬ 
known,  and  in  order  to  use  (4)- (5)  in  practice,  we  also  use  the  following 
approximation 

1  n 

v  «  V  =  —  Sl/fe. ' h  -  Vi] [/ft.  h  -  Vif  ■  (6) 

P  3= 1 

For  the  case  when  the  observation  system  is  scalar,  i.e.,  M  —  1,  the  MxM 
matrix  V  reduces  to  a  scalar  variance  erg,  and  the  equation  (5)  reduces  to 
the  standard  formula 

SEfc  =  \A2C XTx)kk,  k  =  1,2,  ...,p,  (7) 

with  x(9)  an  n  x  p  sensitivity  matrix  for  our  model  given  by 

=  df{tj,9) 
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For  this  scalar  observation  system,  the  approximation  a2  to  cTq  in  (7)  is 
usually  [16]  taken  as 


On 


a  = 


n 


P  3=1 


(9) 


2.1  Traditional  Sensitivity  Functions 

Traditional  sensitivity  functions  (TSF)  are  classical  sensitivity  functions  used 
in  mathematical  modeling  to  investigate  variations  in  the  output  of  a  model 
resulting  from  variations  in  the  parameters  and  the  initial  conditions. 

In  order  to  quantify  the  variation  in  the  state  variable  x(t)  with  respect  to 
changes  in  the  parameter  9  and  the  initial  condition  x(to),  we  are  naturally  led 
to  consider  (traditional)  sensitivity  functions  (TSF)  defined  by  the  derivatives 

Ox 

S0k  (*)  =  a^(*)’  k  =  1,---,p,  (10) 

and 

Ox 

=  l  =  (11) 

Ox  oi 

where  Xo i  is  the  Z-th  component  of  the  initial  condition  xq.  If  the  function  g 
is  sufficiently  regular,  the  solution  x  is  differentiable  with  respect  to  9k  and 
Xq i,  and  therefore  the  sensitivity  functions  sgk  and  rXQl  are  well  defined. 

Often  in  practice,  the  model  under  investigation  is  simple  enough  to  allow 
us  to  combine  the  sensitivity  functions  (10)  and  (11),  as  is  the  case  with 
the  logistic  growth  population  example  discussed  below.  However,  when 
one  deals  with  a  more  complex  model,  as  with  the  agricultural  production 
network  example,  it  is  often  preferable  to  consider  these  sensitivity  functions 
separately  for  clarity  purposes. 

Because  they  are  defined  by  partial  derivatives  which  have  a  local  char¬ 
acter,  the  sensitivity  functions  are  also  local  in  nature.  Thus  sensitivity  and 
insensitivity  (sgk  =  dx/dOk  very  close  to  zero)  depend  on  the  time  interval, 
the  state  values  x ,  and  the  values  of  9  for  which  they  are  considered.  Thus 
for  example  in  a  certain  time  subinterval  we  might  fold  sek  small  so  that  the 
state  variable  x  is  insensitive  to  the  parameter  6k  on  that  particular  inter¬ 
val.  The  same  function  sgk  can  take  large  values  on  a  different  subinterval, 
indicating  to  us  that  the  state  variable  x  is  very  sensitive  to  the  parameter 
6k  on  the  latter  interval.  From  the  sensitivity  analysis  theory  for  dynamical 
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systems,  one  finds  that  s  =  (s^, . . .  ,sq  )  is  an  N  x  p  vector  function  that 
satisfies  the  ODE  system 


s(t)  =  9x(t,x(t),0)s(t)  +go(t,x(t),0), 
s(t0)  0  Nxpi 


(12) 


so  that  the  dependence  of  s  on  (t,  x(t))  as  well  as  9  is  readily  apparent.  Here 
we  denote  by  gx  =  dg/dx  and  by  go  =  dg/dO  the  derivatives  of  g  with  respect 
to  x  and  9 ,  respectively. 

In  a  similar  manner,  the  sensitivity  functions  with  respect  to  the  com¬ 
ponents  of  the  initial  condition  Xo  define  an  N  x  N  vector  function  r  = 
(rx oi,  •  •  • ,  rX0N),  which  satisfies 


(13) 


The  equations  (12)  and  (13)  are  used  in  conjunction  with  equation  (1) 
to  numerically  compute  the  sensitivities  s  and  r  for  general  cases  when  the 
function  g  is  sufficiently  complicated  to  prohibit  a  closed  form  solution  by 
direct  integration. 

Because  the  parameters  may  have  different  units  and  the  state  variables 
may  have  varying  orders  of  magnitude,  sometimes  in  practice  it  is  more 
convenient  to  work  with  the  normalized  version  of  the  TSF,  referred  to  as 
relative  sensitivity  functions  (RSF).  However,  since  in  this  paper  we  are  using 
the  standard  error  approach  to  analyze  the  performance  of  the  least  squares 
algorithm  in  estimating  the  true  parameter  values,  we  will  focus  solely  on 
the  non-scaled  sensitivities,  i.e. ,  TSF. 

2.2  Generalized  Sensitivity  Functions 

Generalized  sensitivity  functions  were  proposed  by  Thomaseth  and  Cobclli 
[26]  as  a  new  tool  in  identification  studies  to  analyze  the  distribution  of  the 
information  content  (with  respect  to  the  model  parameters)  of  the  output 
variables  of  a  system  for  a  given  set  of  observations. 

For  a  scalar  observation  model  with  discrete  time  measurements  (i.e., 
when  M  —  1  and  C  is  a  1  x  N  array  in  (2)),  the  generalized  sensitivity 
functions  (GSF)  are  defined  as 


(14) 
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where  {ti},  l  —  1, ...  ,n  are  the  times  when  the  measurements  are  taken,  and 


1 

F  =  Ej37 (15) 

j=  1  ' 

is  the  corresponding  Fisher  information  matrix.  The  symbol  represents 
element-by-element  vector  multiplication  and,  for  motivation  and  details 
which  led  to  the  definition  above,  the  interested  reader  may  consult  [11,  26]. 
The  Fisher  information  matrix  measures  the  information  content  of  the  data 
corresponding  to  the  model  parameters.  In  (14)  we  see  that  this  information 
is  transferred  to  the  GSF,  making  them  appropriate  tools  to  indicate  the 
relevance  of  the  measurements  in  parameter  estimation  problems. 

We  note  that  the  generalized  sensitivity  functions  (14)  are  vector- valued 
functions  with  the  same  dimension  as  9.  The  fc-th  component  gsk  of  the 
vector  function  gs  represents  the  generalized  sensitivity  function  with  re¬ 
spect  to  9k ■  The  GSF  in  (14)  are  defined  only  at  the  discrete  time  points 
{tj,  j  =  1, . . . ,  n}  and  they  are  cumulative  functions  involving  at  time  ti  only 
the  contributions  of  those  measurements  up  to  and  including  tf,  thus  gsk 
calculates  the  influence  of  measurements  up  to  ti  on  the  parameter  estimate 
for  9k. 

It  is  readily  seen  from  the  definition  that  all  the  components  of  gs  are 
one  at  the  final  time  point  tn,  i.e. ,  gs (tn)  =  1.  If  one  defines  gs(t)  =  0 
for  t  <  ti  (naturally,  gs  is  zero  when  no  measurements  are  collected),  then 
each  component  of  gs  varies  between  0  and  1.  As  developed  in  [26],  the 
time  subinterval  during  which  the  change  in  gsk  has  the  sharpest  increase 
corresponds  to  the  observations  which  provide  the  most  information  in  the 
estimation  of  9k-  That  is,  regions  of  sharp  increases  in  gsk  indicate  a  high 
concentration  of  information  in  the  data  about  9k. 

The  numerical  implementation  of  the  generalized  sensitivity  functions 
(14)  is  straightforward,  since  the  gradient  of  /  with  respect  to  9  (or  xq)  is 
simply  the  Jacobian  of  x  with  respect  to  9  (or  Xo )  multiplied  by  the  obser¬ 
vation  operator  C .  These  Jacobian  matrices  can  be  obtained  by  numerically 
solving  the  sensitivity  ODE  system  (12)  or  (13)  coupled  with  the  system  (1). 
We  will  use  this  approach  to  compute  the  GSF  for  the  agricultural  production 
model.  For  the  other  example,  the  logistic  model,  the  right  side  of  equation 
(1)  is  sufficiently  simple  to  permit  a  closed  form  solution. 


3  The  Verhulst-Pearl  Population  Model 


We  begin  by  considering  the  Verhulst-Pearl  logistic  growth  equation  [22], 
which  approximates  the  evolution  of  a  population  size  over  time  and  is  given 
by 


dx 

dt 


rx 


(16) 


The  constants  K  and  r  represent  the  carrying  capacity  and  the  intrinsic 
growth  rate  respectively.  The  solution  x(t)  of  (16),  representing  the  popula¬ 
tion  number  at  time  t,  is  given  by 


x{t) 


K 


(17) 


where  Xq  =  x(0)  is  the  initial  population  size.  The  solution  x(t)  approaches 
an  asymptote  at  x  =  K  as  t  —>  oo;  this  is  depicted  in  Figure  1. 

The  Verhulst-Pearl  logistic  equation  is  a  relatively  simple  example  with 
easily  studied  dynamics  that  is  useful  in  demonstrating  the  utility  of  the 
traditional  sensitivity  functions  as  well  as  the  generalized  sensitivity  functions 
in  inverse  problems  (see  [9]  for  more  discussions  on  TSF  for  this  system). 
Unless  data  is  sampled  from  regions  with  changing  dynamics,  it  is  possible 
that  some  of  the  parameters  will  be  difficult  to  estimate.  Moreover,  the 
parameters  that  are  obtainable  may  have  high  standard  errors  as  a  result 
of  introducing  redundancy  in  the  sampling  region.  In  order  to  demonstrate 
this  for  the  logistic  growth  problem,  we  will  examine  varying  behavior  in 
the  model  depending  on  the  region  from  which  tj  is  sampled.  We  consider 
points  T\  and  r2,  as  depicted  in  Figure  1,  partitioning  the  curve  into  three 
distinct  regions:  0  <  tj  <  Ti,  T\  <  tj  <  r2,  and  r2  <  tj  <  T,  with  T 
sufficiently  large  for  our  solution  to  be  near  its  asymptote  x  =  K.  Based  on 
the  changing  dynamics  of  the  curve  in  Figure  1,  we  expect  differences  in  the 
ability  to  estimate  parameters  depending  on  the  region  in  which  the  solution 
is  observed. 

To  illustrate  these  ideas,  we  carried  out  estimation  procedures  for  the 
parameters  6  =  (K,  r,  x0)  in  the  logistic  growth  population  model  using  or¬ 
dinary  least  square  procedures  with  both  numerically  generated  (no-noise) 
and  noisy  simulated  “data”.  We  produced  data  sets  { Vy }y= i  by  evaluating 
the  numerical  solution  of  (16)  with  the  “true”  value  parameters  6q  at  tj  and 
then  adding  random  noise  in  some  cases.  With  a  small  dimension  parameter 
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K 


Figure  1:  Distinct  regions  of  growth  in  the  Verhulst-Pearl  solution  curve. 


space  as  in  this  example,  a  Nelder-Mead  optimization  algorithm  is  sufficient 
for  the  inverse  problem.  We  hence  use  the  MATLAB  function  fminsearch 
to  minimize  the  cost  functional 


n 


(18) 


with  respect  to  9  and  using  an  initial  guess  9°  for  the  optimization  algorithm. 
In  order  to  avoid  what  is  typically  called  an  inverse  crime,  we  evaluate  /  in 
the  cost  function  using  the  ODE  solver  odel5s  with  (16),  which  returns  the 
numerical  approximation  to  the  solution  fit,  9)  =  x(t;  K,r,x0),  rather  than 
using  the  analytical  solution  (17).  This,  in  effect,  produces  “noisy  data”  even 
when  no  additional  noise  is  added.  In  this  note  we  illustrate  ideas  with  a 
specific  example,  taking  90  =  (K,r,xo)  =  (17.5,0.7,0.1);  the  corresponding 
solution  curve  x(t)  can  be  seen  in  Figure  2(a). 

3.1  Traditional  Sensitivities 

We  consider  an  ordinary  least  squares  problem  for  the  estimation  of  the 
parameters  9  =  ( K ,  r,  Xo)  in  the  logistic  growth  model  (16),  using  the  explicit 
solution  given  by  (17)  and  then  examining  the  sensitivities  with  respect  to 
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the  parameters  K,  r,  and  x$.  We  can  readily  compute  the  partial  derivatives 


dx  Xq(1  —  e~rt) 

dK  (x0  +  (K  —  x0)e~rt)2  ’ 
dx  Kx0(K  —  x0)te~rt 

dr  (xq  +  (. K  —  xo)e~rt )2  ’ 
dx  K2e~rt 

dx o  (x0  + (K  -  x0)e-rt)2' 


which  are  the  traditional  sensitivity  functions  sk,  sr  and  sxo. 

We  analyze  the  TSF  corresponding  to  each  parameter  in  the  initial  region 
of  the  curve,  where  the  solution  approaches  xq.  When  we  consider  the  initial 
region  of  the  curve,  where  0  <  tj  <  T\  for  j  —  1, . . . ,  n,  we  have 


dx(tj ) 
dK 


0, 


dx(tj ) 
dr 


0, 


dx(tj ) 
dx  o 


1; 


this  follows  from  considering  the  limits  of  the  sensitivity  functions  (19)  as 
t  — >  0.  Based  on  the  above  analytical  findings,  which  indicate  low  sensitivities 
with  respect  to  K  and  r,  we  expect  to  have  little  ability  to  determine  these 
parameters  when  we  sample  data  from  [0,7i];  however  we  should  be  able  to 
estimate  ^o- 

We  next  consider  the  region  of  the  curve  which  is  near  the  asymptote  at 
x  =  K ,  in  this  case  for  r2  <  tj  <  T,  j  —  1, . . . ,  n.  Here  we  find  that  by 
considering  the  limits  as  t  — >  oo,  we  have  the  approximations 


dx(tj ) 
dK 


dx(tj ) 
dr 


0, 


dx(tj ) 
dx  o 


0. 


Based  on  these  approximations,  we  expect  to  be  able  to  estimate  K  well  when 
we  sample  data  from  [r2,  T] .  However,  using  data  only  from  this  region,  we 
do  not  expect  to  be  able  to  estimate  Xq  or  r. 

Finally,  we  consider  the  part  of  the  solution  curve  where  T\  <  tj  <  r2  for 
j  =  1 , ,n  and  where  it  has  nontrivially  changing  dynamics.  We  note  that 
the  partial  derivative  values  differ  greatly  from  the  values  in  regions  [0,  Ti] 
and  [t2,T],  When  [ti,t2]  is  included  in  the  sampling  region  we  expect  to 
recover  good  estimates  for  all  three  parameters. 

Our  analytical  observations  are  fully  consistent  with  information  con¬ 
tained  in  the  graphs  of  the  TSF  illustrated  in  Figure  2(b)  for  T  =  25.  We 
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x(t) 


Model  Sensitivities 


Figure  2:  (a)  Logistic  curve  (b)  TSF  corresponding  to  each  parameter  for 
the  logistic  curve  with  6  =  (17.5,0.7,0.1). 


note  that  the  curve  Sk  slowly  increases  with  time  and  it  appears  that  the 
solution  is  insensitive  to  K  until  around  the  flex  point  of  the  logistic  curve, 
which  occurs  shortly  after  t  —  7  in  this  case.  The  sensitivities  sk  and  sr  both 
are  close  to  zero  when  t  is  near  the  origin,  and  hence  we  deduce  that  both 
K  and  r  will  be  difficult  or  impossible  to  obtain  using  data  in  that  region. 
Also,  we  observe  that  sXo  and  sr  are  nearly  zero  in  [15,25],  which  suggests 
that  we  will  be  unable  to  estimate  xo  or  r  using  observations  in  that  region. 

In  order  to  computationally  illustrate  how  the  traditional  sensitivity  the¬ 
ory  applies  to  our  logistic  growth  example,  we  present  here  a  summary  of 
findings  obtained  using  numerically  (no-noise  added)  simulated  data  gener¬ 
ated  with  odel5s.  We  also  restricted  the  time  domain  from  which  we  sampled 
data  points  to  one  of  the  three  regions  of  interest  [0,  tl],  [ti,  72]  or  [72,  T]  for 
each  inverse  problem  calculation.  Then  we  determined  the  standard  errors 
of  each  parameter  set  in  order  to  analyze  the  success  of  the  algorithm  at  ap¬ 
proximating  the  various  components  of  6.  (Other  computations  with  varying 
levels  of  added  noise  in  the  “data”  are  presented  in  [9]  where  findings  are 
similar  to  those  reported  here.) 

We  first  sampled  data  from  the  region  where  the  solution  curve  is  close 
to  a’o,  and  for  this  example  we  considered  the  region  [0, 1].  We  used  several 
initial  guesses  in  our  MATLAB  solver,  and  expected  with  each  guess  that  we 
would  be  able  to  achieve  an  appropriate  estimate  for  x0j  but  perhaps  not  for 
r  or  K .  In  actuality  we  were  able  to  obtain  close  estimates  for  both  r  and 
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Xq,  as  reported  in  Table  1. 


0° 

6 

Standard  Errors 

Values  obtained  using  data  in  0  <  t  <  1 

(8, 1,0.3) 
(15,2,0.2) 
(30,0.3,0.5) 
(10,0.1,0.3) 

(10.686,0.7038,0.09998) 

(26.8241,0.6978,0.1000) 

(34.1993,0.6970,0.1000) 

(16.9140,0.7001,0.1000) 

(10.398, 0.003589, 3.4206e-05) 
(7.4758, 0.002581, 2.4593e-05) 
(9.9947, 0.003450, 3.2879e-05) 
(0.7924, 0.000274, 2.6068e-06) 

Values  obtained  using  data  in  1  <  t  <  15 

(8, 1,0.3) 
(15,2,0.2) 
(30,0.3,0.5) 
(10,0.1,0.3) 

(17.4985,0.6998,0.1001) 

(17.4985,0.6998,0.1001) 

(17.4985,0.6998,0.1001) 

(17.4984,0.6998,0.1001) 

(0.002062,0.0003049,0.0002142) 

(0.002063,0.0003050,0.0002143) 

(0.002061,0.0003048,0.0002141) 

(0.002063,0.0003051,0.0002144) 

Values  obtained  using  data  in  15  <  t  <  25 

(8, 1,0.3) 
(15,2,0.2) 
(30,0.3,0.5) 
(10,0.1,0.3) 

(17.4879,1.1427,0.0495) 
(17.4877,1.7838,0.1908) 
(17.5017,0.5887,0.5374) 
(17.5023,0.5468,1 .0098) 

(0.02912,1.5016,2.3000) 

(0.03076,1.5865,2.4301) 

(0.00268,0.1381,0.2115) 

(0.00328,0.1689,0.2587) 

Table  1:  The  optimized  8  values  and  corresponding  standard  errors,  on  the 
given  intervals  with  60  =  (17.5,  0.7,  0.1)  implemented  using  a  computationally 
noisy  data  set. 


Upon  further  examination  we  found  that  r  and  Xq  are  highly  correlated 
when  0  <  t  <  1,  which  is  evident  by  the  magnitude  of  their  correlation 
coefficients,  given  in  Table  2.  By  studying  each  iteration  of  the  MATLAB 
solver,  we  observed  that  a  good  estimate  for  xq  was  easily  obtained,  and  then, 
due  to  the  high  correlation  between  r  and  Xq,  r  was  eventually  obtained  each 
time.  However,  true  to  our  predictions,  K  was  never  reasonably  estimated 
using  data  from  this  region. 

Recall  that  correlation  coefficients  range  from  —1  to  1,  where  higher  mag¬ 
nitudes  indicate  stronger  correlation  between  the  corresponding  parameters. 
The  correlation  coefficients  displayed  in  Table  2  were  computed  using  the 
standard  definition 

cov(X: 

CTx&Y 

where  X  and  Y  are  two  random  variables  with  mean  fix  and  fly  and  variances 


corr(X,  Y ) 
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K 

r 

x0 

K 

1.0000 

-0.3434 

0.2855 

r 

-0.3434 

1.0000 

-0.977 

x0 

0.2855 

-0.977 

1.0000 

Table  2:  Correlation  coefficients  for  parameters 


o\  and  cry.  The  parameter  estimates  6n  can  be  interpreted  as  realizations  of 
a  random  variable  0n,  with  distribution  which,  from  the  asymptotic  theory 
of  statistical  analysis  for  the  least  square  algorithm,  can  be  approximated  by 
a  normal  distribution  as  n  — *>  oo,  i.e. ,  for  n  large,  071  ■A/'p(6'0,crg(xTx)  x)  is 

a  good  approximation.  The  correlation  coefficients  between  two  components 
of  9  are  thus  given  by 

cov(9k,  6 1 
(?ek°el 

Using  the  definition  of  the  covariance  matrix,  we  have  that  cov{9k ,  Of)  is 
simply  the  (k,l)-th  element  of  a q(xTx)~1  and  the  standard  deviations  agk 
and  agl  are  the  square  roots  of  the  (k,  k)- th  and  (l,  l)- th  diagonal  entries. 
This  was  used  to  compute  the  approximate  correlation  coefficients  of  Table  2. 

The  logistic  model  above  was  considered  and  analyzed  in  [9]  using  tra¬ 
ditional  sensitivity  functions.  However  the  model  was  formulated  with  a 
different  parameterization:  instead  of  9  —  (K,r,x o),  the  parameter  set  was 
given  by  (3  —  (a,  b,  x0)  =  (r,  j^,x0).  Studying  the  TSF  curves  with  this  other 
parameterization,  we  found  there  was  no  difficulty  in  predicting  the  regions  in 
which  the  state  was  sensitive  to  each  parameter  with  no  consideration  given 
to  the  correlation  coefficients.  With  a  new  parameterization,  6  =  (K,r,x o), 
the  correlation  between  the  coefficients  reveals  direct  information  on  our  abil¬ 
ity  to  predict  the  identihability  for  each  parameter.  Comparing  these  findings 
with  those  of  [9] ,  we  see  that  the  parameterization  of  the  model  is  critical  in 
sensitivity  studies. 

We  next  considered  the  region  where  the  dynamics  of  the  curve  are  chang¬ 
ing,  and  for  our  example  this  region  is  the  interval  [1, 15].  As  expected,  re¬ 
gardless  of  the  initial  guess  that  is  used  in  the  solver,  we  can  generally  obtain 
reasonable  estimates  for  K,  r  and  Xq]  this  can  be  seen  in  Table  1. 

Finally  we  sampled  data  from  the  region  where  the  curve  approaches  an 
asymptote  at  K,  and  for  this  example  we  considered  the  region  [15,25].  We 


corr(Ok,  9t ) 
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used  several  initial  guesses  in  our  MATLAB  solver,  and  expected  that  with 
each  guess  we  would  be  able  to  achieve  an  appropriate  estimate  for  K,  but 
not  for  r  or  xq.  This  is  precisely  what  occurred  with  every  initial  guess. 

We  see  that  the  TSF,  used  in  conjunction  with  the  correlation  coefficients, 
provided  sufficient  information  in  these  examples  to  predict  which  parame¬ 
ters  could  be  determined  using  data  from  different  regions.  Similar  results 
obtained  using  noise-added  data  with  the  logistic  model,  along  with  some 
pitfalls  of  the  TSF,  can  be  found  in  [9]. 

3.2  Generalized  Sensitivities 

We  next  illustrate  the  utility  of  the  generalized  sensitivity  functions  by  ap¬ 
plying  the  theory  to  the  logistic  growth  model  (16).  We  start  by  numerically 
computing  the  GSF  using  equation  (14)  with  a  =  1  and  the  true  value  pa¬ 
rameters  60  =  (17.5,0.7,0.1).  The  plots  of  these  functions  are  shown  in 
Figure  3(b)  where  one  can  observe  obvious  regions  of  steep  increase  in  each 
curve.  For  the  curves  gsXQ(t),  gsr(t )  and  gsK(t),  we  find  by  visual  inspection 

x(t)  GSF 


Figure  3:  (a)  Logistic  curve  (b)  GSF  corresponding  to  each  parameter  for 
the  logistic  curve  with  6  =  (17.5,0.7,0.1). 


that  these  regions  are  approximately  [4.5,  7.5],  [7,11]  and  [12,25],  respec¬ 
tively.  By  the  aforementioned  generalized  sensitivity  theory,  if  we  increase 
the  number  of  data  points  sampled  in  one  of  these  regions,  the  estimation  of 
the  corresponding  parameter  is  expected  to  improve. 

In  order  to  implement  the  GSF  for  the  logistic  growth  example,  we  intro¬ 
duced  noise  into  the  simulated  data  by  adding  Gaussian  noise  0,  Uq) 
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for  j  =  1  To  obtain  the  results  in  Tables  3,  4,  and  5,  we  used  a 

randomly  generated  noisy  data  set  with  <To  =  0.5.  We  initially  sampled  n 


n 

rn 

t-GSF 

tnonGSF 

50 

0 

(0.10748,0.029371,0.021604) 

(0.10748,0.029371,0.021604) 

50 

25 

(0.077068,0.027991,0.020752) 

(0.10441,0.021464,0.015622) 

50 

50 

(0.059962,0.025752,0.019184) 

(0.099866,0.01754,0.012632) 

125 

0 

(0.061814,0.016807,0.01236) 

(0.061814,0.016807,0.01236) 

125 

25 

(0.053069,0.016569,0.012234) 

(0.061051,0.014428,0.010561) 

125 

50 

(0.046431,0.016095,0.011919) 

(0.06453,0.013726,0.010003) 

125 

125 

(0.037958,0.016276,0.012124) 

(0.061772,0.010767,0.0077615) 

250 

0 

(0.047175,0.012805,0.0094158) 

(0.047175,0.012805,0.0094158) 

250 

50 

(0.039627,0.012356,0.009123) 

(0.046264,0.010884,0.0079648) 

250 

100 

(0.035866,0.012422,0.009198) 

(0.045259,0.0095868,0.0069842) 

250 

250 

(0.027388,0.011737,0.0087426) 

(0.04239,0.0073602,0.0053052) 

n 

rn 

t uniform 

50 

0 

(0.10748,0.029371,0.021604) 

50 

25 

(0.086565,0.023623,0.017374) 

50 

50 

(0.077352,0.021049,0.01548) 

125 

0 

(0.061814,0.016807,0.01236) 

125 

25 

(0.057056,0.015504,0.011401) 

125 

50 

(0.053388,0.014502,0.010664) 

125 

125 

(0.047175,0.012805,0.0094158) 

250 

0 

(0.047175,0.012805,0.0094158) 

250 

50 

(0.043685,0.011854,0.0087166) 

250 

100 

(0.0404,0.01096,0.0080592) 

250 

250 

(0.034008,0.0092229,0.0067818) 

Table  3:  This  table  shows  standard  errors  corresponding  to  the  addition  of 
m  extra  points  in  or  excluding  [12,25],  as  it  is  the  period  of  steepest  increase 
for  the  parameter  K. 


data  points  uniformly  over  the  entire  region  [0,25],  and  then  we  sampled  an 
additional  m  points  from  varying  parts  of  the  entire  region,  comparing  the 
standard  errors  of  each  trial. 

Each  parameter  in  9  has  a  specific  interval  of  steepest  increase  according 
to  the  corresponding  GSF,  and  the  m  additional  points  are  sampled  according 
to  those  regions.  In  each  of  Tables  3,  4, and  5,  we  consider  three  separate  time 
grids: 
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•  tcsF  consists  of  n  uniform  time  points  over  [0,25]  with  m  points  added 
in  the  area  of  steepest  increase  according  to  the  GSF  for  each  parame¬ 
ter, 

•  tnonGSF  consists  of  71  uniform  time  points  over  [0,25]  with  m  points 
added  everywhere  except  the  area  of  steepest  increase  according  to  the 
GSF  for  each  parameter, 

•  tuniform.  is  the  time  grid  referring  to  n  uniform  time  points  over  [0,25] 
with  m  additional  points  added  as  uniformly  as  possible  over  the  entire 
region. 

For  example,  in  Table  3,  the  tcsF  column  refers  to  the  standard  errors 
generated  by  optimizing  9  over  n  uniformly  distributed  data  points  in  [0,25] 
with  m  additional  data  points  sampled  from  the  steepest  increase  interval 
corresponding  to  K:  [12,25].  The  standard  errors  in  the  tnonGSF  column 
refer  to  the  same  n  initial  data  points,  with  the  additional  m  points  sampled 
from  outside  the  GSF  suggested  region  for  K. 

In  general,  values  of  standard  errors  are  meaningful  only  relative  to  the 
estimated  values  of  the  corresponding  parameters.  Thus,  one  should  report 
both  the  estimated  parameter  values  and  the  corresponding  SE.  However, 
here  and  in  Section  4,  we  report  only  the  changes  in  SE  as  data  sets  are 
changed.  This  is  because  our  primary  focus  is  on  these  changes  in  SE  and 
because  the  corresponding  OLS  estimates  are  near  (same  order  magnitude) 
the  true  values  90;  hence  it  suffices  to  simply  report  only  90  and  the  changes 
in  SE. 

Note  that  since  9  =  (K,  r,  x0),  the  standard  error  for  K  is  indicated  as  the 
first  entry  in  each  of  the  ordered  sets  in  each  table,  i.e.,  SEk  =  SEg1,  and 
Tabic  3  refers  to  the  steepest  region  of  increase  in  the  GSF  corresponding  to 
K .  Similarly  Tables  4  and  5  use  the  GSF  regions  corresponding  to  r  and  xq, 
respectively,  and  so  SEq2  and  then  SEg 3  are  the  entries  of  interest. 

We  would  have  expected  that  for  each  parameter,  the  corresponding  stan¬ 
dard  error  in  the  tasF  column  would  generally  improve  as  m  increased.  In 
actuality,  the  corresponding  standard  errors  generally  improved  when  m  ad¬ 
ditional  points  were  sampled  in  all  three  grids  for  each  parameter.  However, 
we  get  better  results  with  the  tcsF  grid  than  with  the  tuniform  grid,  and  the 
tnonGSF  grid  was  the  worst  of  the  three.  Hence  we  see  that  by  catering  to  the 
GSF-reconnnended  regions,  we  are  able  to  obtain  better  standard  errors  for 
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n 

rn 

tGSF 

t-nonGSF 

50 

0 

(0.10748,0.029371,0.021604) 

(0.10748,0.029371,0.021604) 

50 

25 

(0.097735,0.021034,0.016379) 

(0.085987,0.025298,0.01814) 

50 

50 

(0.10215,0.019931,0.015575) 

(0.071516,0.022395,0.015627) 

125 

0 

(0.061814,0.016807,0.01236) 

(0.061814,0.016807,0.01236) 

125 

25 

(0.061485,0.014585,0.011165) 

(0.055069,0.0155,0.011258) 

125 

50 

(0.063405,0.013979,0.010842) 

(0.052744,0.015359,0.011012) 

125 

125 

(0.062507,0.012136,0.009478) 

(0.043153,0.013462,0.0093937) 

250 

0 

(0.047175,0.012805,0.0094158) 

(0.047175,0.012805,0.0094158) 

250 

50 

(0.046673,0.011057,0.0084657) 

(0.043127,0.012158,0.0088167) 

250 

100 

(0.045837,0.010087,0.0078221) 

(0.038498,0.01117,0.0080133) 

250 

250 

(0.043884,0.0085094,0.0066447) 

(0.032496,0.010137,0.0070706) 

n 

rn 

^ uniform 

50 

0 

(0.107480,0.029371,0.021604) 

50 

25 

(0.086565,0.023623,0.017374) 

50 

50 

(0.077352,0.021049,0.015480) 

125 

0 

(0.061814,0.016807,0.012360) 

125 

25 

(0.057056,0.015504,0.011401) 

125 

50 

(0.053388,0.014502,0.010664) 

125 

125 

(0.047175,0.012805,0.009416) 

250 

0 

(0.047175,0.012805,0.0094158) 

250 

50 

(0.043685,0.011854,0.0087166) 

250 

100 

(0.040400,0.010960,0.0080592) 

250 

250 

(0.034008,0.009223,0.0067818) 

Table  4:  This  table  shows  standard  errors  corresponding  to  the  addition  of 
m  extra  points  in  or  excluding  [7,11],  as  it  is  the  period  of  steepest  increase 
for  the  parameter  r. 


the  corresponding  parameters  than  if  we  had  merely  increased  the  number 
of  points  sampled  over  the  entire  region. 

3.3  “Forced-to-one”  Artifact 

We  first  note  that  the  shape  of  the  TSF  curves  remains  the  same  regardless 
of  the  amount  of  data  that  is  sampled,  whereas  the  GSF  curves  are  data 
dependent  and  change  shape  with  varying  amounts  of  data.  We  can  see  in 
Figure  4  that  when  we  restrict  the  data  set  for  the  logistic  growth  model  to 
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n 

rn 

t-GSF 

t-nonGSF 

50 

0 

(0.10748,0.029371,0.021604) 

(0.10748,0.029371,0.021604) 

50 

25 

(0.10338,0.021906,0.014617) 

(0.081923,0.02323,0.017621) 

50 

50 

(0.099638,0.018577,0.012139) 

(0.07291,0.021106,0.016235) 

125 

0 

(0.061814,0.016807,0.01236) 

(0.061814,0.016807,0.01236) 

125 

25 

(0.064238,0.015175,0.010452) 

(0.056836,0.015457,0.011486) 

125 

50 

(0.059416,0.012951,0.0086999) 

(0.052824,0.01476,0.011091) 

125 

125 

(0.05867,0.010902,0.0071197) 

(0.043777,0.012942,0.0099389) 

250 

0 

(0.047175,0.012805,0.0094158) 

(0.047175,0.012805,0.0094158) 

250 

50 

(0.046931,0.011043,0.0076218) 

(0.043069,0.011892,0.0088565) 

250 

100 

(0.046695,0.010136,0.006819) 

(0.040509,0.011467,0.0086329) 

250 

250 

(0.043801,0.0081056,0.0052964) 

(0.031893,0.0095317,0.0073385) 

n 

rn 

t uniform 

50 

0 

(0.10748,0.029371,0.021604) 

50 

25 

(0.086565,0.023623,0.017374) 

50 

50 

(0.077352,0.021049,0.01548) 

125 

0 

(0.061814,0.016807,0.01236) 

125 

25 

(0.057056,0.015504,0.011401) 

125 

50 

(0.053388,0.014502,0.010664) 

125 

125 

(0.047175,0.012805,0.0094158) 

250 

0 

(0.047175,0.012805,0.0094158) 

250 

50 

(0.043685,0.011854,0.0087166) 

250 

100 

(0.04040,0.010960,0.0080592) 

250 

250 

(0.034008,0.0092229,0.0067818) 

Table  5:  This  table  shows  standard  errors  corresponding  to  the  addition  of  m 
extra  points  in  or  excluding  [4. 5, 7.5],  as  it  is  the  period  of  steepest  increase 
for  the  parameter  xq- 


[0,  2],  the  TSF  curves  look  the  same  as  those  where  we  merely  zoom  in  to  [0, 2] 
after  sampling  data  over  the  entire  region  [0,  25].  Then  in  Figure  5  the  GSF 
curves  look  very  different  when  the  sampled  data  is  only  from  the  interval 
[0,  2]  as  compared  to  when  the  data  is  sampled  from  [0,  25]  and  zoomed- in 
to  [0,2],  If  an  insufficient  amount  of  data  is  used  for  parameter  estimations, 
the  GSF  curves  can  be  misleading  because  by  definition  the  GSF  are  forced 
to  be  equal  to  one  by  the  end  of  the  data  set.  This  so-called  “forced-to-one” 
artifact  can  cause  misleading  regions  of  steep  increase  in  the  GSF  curves  as 
can  be  seen  in  the  curve  gsx  when  t  E  [1.7, 2]  in  Figure  5(a).  Observe  that  in 


19 


TSF  on  [0,2] 


TSF  on  [0,25],  zoomed  in  to  [0,2] 


Figure  4:  In  (a)  and  (b)  we  see  the  TSF  of  the  parameters  when  sampling 
data  in  [0,2]  and  the  zoomed-in  portion  when  data  was  observed  from  the 
entire  region  [0,25]. 


GSF  on  [0,2] 


GSF  on  [0,25],  zoomed  in  to  [0,2] 


Figure  5:  In  (a)  we  see  the  GSF  of  the  parameters  when  sampling  data  in 
[0,2]  and  in  (b)  we  have  the  zoomed-in  portion  [0,2]  when  data  was  observed 
from  the  entire  region  [0,25]. 


Figure  4  we  can  see  that  the  state  is  clearly  not  sensitive  to  the  K  parameter 
in  region  [0,2],  and  hence  sampling  additional  data  points  in  the  period  of 
misleading  increase  would  merely  make  estimates  worse,  with  an  increase  in 
standard  error.  We  can  also  see  another  example  of  this  in  Figure  6  where 
we  compare  the  GSF  curves  generated  from  data  sampled  from  [0, 0.2]  to  the 
curves  generated  when  data  is  sampled  from  [0,  25].  In  Figure  6(b)  it  is  clear 
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that  there  is  no  significant  period  of  increase  for  either  K  or  r,  however  in 
(a)  it  appears  that  the  corresponding  GSF  curves  both  have  distinct  regions 
of  (again  in  this  case  misleading)  increase.  Therefore,  it  is  important  to  note 
that  while  the  GSF  curves  can  be  misleading  when  a  limited  portion  of  data  is 
obtainable,  the  TSF  curves  can  still  be  used  appropriately  in  such  sensitivity 
studies. 


Figure  6:  In  (a)  we  see  the  GSF  of  the  parameters  when  sampling  data 
in  [0,0.2]  and  in  (b)  we  have  the  zoomed-in  portion  [0,0.2]  when  data  was 
observed  from  the  entire  region  [0,25]. 


It  is  also  important  to  remark  on  the  Fisher  information  matrix  which 
appears  in  the  definition  of  generalized  sensitivity  functions  (14).  The  in¬ 
formation  on  the  parameters  provided  by  the  measurement  ?/*  is  quantified 
by  the  derivative  of  an  information  index  with  respect  to  When  we  try 
to  estimate  the  parameters  r  and  Xq  with  data  from  the  interval  [15,  25] 
alone  (see  Table  1),  we  obtain  large  errors  which  increase  in  magnitude  when 
the  number  of  sample  points  increases.  Although  initially  unexpected,  this 
phenomenon  can  be  explained  based  on  the  GSF  discussion  presented  above. 
When  we  sample  additional  data  points  from  the  region  where  the  traditional 
sensitivity  curves  are  flat  and  the  generalized  sensitivity  functions  exhibit  the 
“forced-to-one”  artifact,  we  actually  introduce  redundancy  in  the  sensitivity 
matrix.  This  considerably  increases  the  condition  number  of  the  Fisher  in¬ 
formation  matrix,  which  in  turn,  by  the  Crammer-Rao  inequality,  causes  the 
variance  of  the  unbiased  estimator  to  be  very  large,  making  our  estimates 
less  useful. 
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Although  there  is  improvement  when  the  GSF-recommended  regions  are 
considered,  the  amount  of  additional  points  sampled  to  garner  the  improved 
standard  errors  needs  to  be  taken  into  consideration.  Depending  on  the 
problem,  the  cost  may  be  too  high  to  sample  additional  data  points  for  the 
slightly  improved  results.  However,  in  other  cases  the  additional  sampling 
may  be  worth  the  improvement.  While  a  useful  tool,  the  GSF  may  not  be 
an  efficient  choice  in  some  parameter  improvement  attempts. 

4  An  Agricultural  Production  Model 

We  continue  our  investigations  on  the  relevance  of  the  generalized  sensitiv¬ 
ity  functions  to  parameter  estimations  problems,  in  this  case  with  a  more 
complex  example  of  an  agricultural  production  network  model  formulated  in 
terms  of  the  nonlinear  dynamical  system 

ci(t)  =  -«iCi(t)(Z2  -  c2(t))  +  +  /c4min(c4(t),  sm) 

c2(t)  =  —  K2c2(t)(h  ~  c3(Z))+  +  KiCi(t)(Z2  -  c2(t))+  (20) 

c3(f)  =  -K3c3(i)(Z4  -  c4(t))+  +  K2c2(t)(/3  -  c3(i))+ 

c4(t)  =  — /c4min (c4(t),  sm)  +  K3c3(t)(l4  -  c4(t))+ 

together  with  the  initial  conditions 

c(0)  =  c0.  (21) 

The  system  (20)  is  the  deterministic  limit  for  large  populations  (in  a 
sense  made  precise  in  [2])  of  a  continuous  time  discrete  state  Markov  Chain 
proposed  in  [2]  to  model  the  flow  and  the  impact  of  eventual  disturbances 
in  a  swine  production  network.  For  simplicity,  the  modeled  network  is  as¬ 
sumed  to  consist  of  four  levels  of  production  nodes:  Growers  (TV) ),  Nurseries 
(N 2),  Finishers  (N3),  and  Processing  Plants  (iV4),  and  each  node  represents 
an  aggregation  of  all  the  production  units  corresponding  to  that  level  in  the 
production  process.  Although  unrealistic,  it  is  assumed  that  there  are  no 
losses  in  the  first  three  nodes  of  the  production  network  and  the  only  deaths 
occur  at  the  processing  plants.  Another  important  assumption  is  that  the 
network  is  closed,  i.e.,  there  is  a  direct  flow  from  node  N4  to  node  N\ ,  instantly 
replenishing  the  network  and  keeping  the  total  population  size  constant  in 
time.  We  note  that  this  assumption  is  realistic  when  the  network  is  efficient 
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and  operates  at  or  near  full  capacity  (i.e.,  when  the  number  of  animals  re¬ 
moved  from  the  chain  are  immediately  replaced  by  new  production/growth, 
avoiding  significant  idle  times). 

The  state  variables  Cj(t),  %  —  1, . . . ,  4,  represent  the  swine  population  size 
(roughly  speaking,  an  ensemble  average  in  the  sense  of  the  Markov  Chain 
model  of  [2])  at  the  nodes  Ni  at  time  t.  The  parameters  «$,  i  =  1, ...  ,4, 
in  (20)  represent  the  transition  rates  between  consecutive  nodes  and  h,  i  — 
2,3,4,  represent  the  maximum  capacitates  at  the  nodes  AC  There  is  no 
capacity  constraint  at  node  N\,  but  there  is  a  maximum  slaughtering  capacity 
at  node  N, 4  (processing  plants)  which  we  denote  by  sm.  For  any  real  z, 
the  symbol  (z)+  is  defined  as  the  positive  part  of  z,  i.e.,  (z)+  =  max(z,0). 
The  numerical  values  of  all  the  parameters  which  we  used  in  our  analysis 
presented  here  are  listed  in  Table  6.  For  more  details  about  the  model  (20) 
and  its  derivation,  the  interested  reader  may  consult  [2] . 


Parameters 

Definition 

Values 

Units 

Ki 

scaled  rate  at  node  1 

2.879 

1 / days 

k2 

scaled  rate  at  node  2 

1.093 

1 / days 

scaled  rate  at  node  3 

1.51 

1 / days 

K  4 

scaled  rate  at  node  4 

1 

1 / days 

k 

scaled  capacity  at  node  2 

2.387-  10"1 

dimensionless 

h 

scaled  capacity  at  node  3 

6.498  •  10"1 

dimensionless 

u 

scaled  capacity  at  node  4 

5.67-  10~2 

dimensionless 

$171 

scaled  slaughter  capacity 

1.039  •  10"1 

dimensionless 

ci(0) 

scaled  initial  condition 

9.45  •  10~2 

dimensionless 

c2(0) 

scaled  initial  condition 

2.221  •  10"1 

dimensionless 

cs(0) 

scaled  initial  condition 

6.313  •  10"1 

dimensionless 

c4(0) 

scaled  initial  condition 

5.19  ■  10"2 

dimensionless 

Table  6:  Aggregated  agricultural  network  model:  Parameters  for  determin¬ 
istic  simulations  (numbers  of  pigs  are  in  thousands). 


We  considered  a  series  of  least  square  problems  using  simulated  noisy  data 
in  order  to  probe  the  utility  of  the  GSF  pertaining  to  estimation  problems. 
Since  we  have  three  categories  of  parameters  in  our  model  (transition  rates, 
capacities  and  initial  conditions  of  the  network)  we  tried  to  estimate  the  pa¬ 
rameters  in  each  category  when  all  the  others  remain  fixed.  The  simulated 
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data  was  generated  by  first  numerically  solving  the  system  (20)  with  param¬ 
eter  and  initial  condition  values  given  in  Table  6  (we  will  refer  to  these  values 
as  the  true  parameter  values  in  these  discussions),  and  then  adding  Gaussian 
noise  of  zero  mean  and  standard  deviation  0.1  to  the  solution  obtained  at 
each  observation  point. 

We  begin  with  the  problem  of  estimating  the  transition  rates  Ki,  k2, 
Ks  and  K4  when  the  l  and  Cok  parameters  are  fixed  at  the  values  given  in 
Table  6.  For  the  true  parameter  values  9q  =  R  =  (2.88,1.09,1.51,1),  we 
plot  the  generalized  sensitivity  functions  (14)  and  the  traditional  sensitivity 
functions  in  Figure  7.  In  both  of  these  plots,  one  can  observe  a  distinct  time 

Generalized  Sensitivity  Functions  with  respect  to  k2,k3,k4  Traditional  Sensitivity  Functions  with  respect  to  k2,k3,k4 


Figure  7:  Generalized  and  Traditional  Sensitivity  Functions  with  respect  to 
Ki,  K2,  k3,  fiq  corresponding  to  parameters  values  given  in  Table  6. 


point  near  t  =  30,  where  the  dynamics  of  both  the  GSF  and  TSF  curves 
change.  Between  0  and  this  point,  the  TSF  plots  exhibit  a  sharp  monotonic 
increase,  and  then  they  reach  a  steady  state  very  quickly.  On  the  contrary, 
the  GSF  curves  increase/decrease  steeply  before  reaching  this  time  point, 
and  then  are  forced  to  one.  According  to  the  GSF  theory,  the  approximate 
interval  [0,  30]  is  the  region  in  which  measurements  are  the  most  informative 
for  estimating  the  true  parameters  R.  This  means  that  by  sampling  additional 
data  points  here,  we  expect  to  obtain  more  information  about  R,  resulting  in 
more  accurate  estimates  for  these  parameters. 

By  comparing  the  GSF  plots,  we  also  observe  that  strong  correlations 
exist  between  «q,  /c2,  k3  and  k4.  The  correlation  coefficients  between  these 
parameters,  given  in  Table  7,  reflect  the  dynamics  of  the  curves  shown  in 
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re  1 

re2 

re3 

K4 

Ki 

1.000 

0.922 

0.879 

0.902 

re  2 

0.922 

1.000 

0.933 

0.910 

re  3 

0.879 

0.933 

1.000 

0.730 

k4 

0.902 

0.910 

0.730 

1.000 

Table  7:  Correlation  coefficients  for  each  parameter 


Data  points 

in 

Standard  Errors  for 

Data  Set 

[0,30] 

[40,210] 

Total 

re  1 

re2 

re3 

k4 

KDS ! 

10 

35 

45 

0.315 

0.108 

0.166 

0.101 

kds2 

30 

35 

65 

0.205 

0.067 

0.114 

0.062 

kds3 

10 

57 

67 

0.244 

0.089 

0.136 

0.085 

KDSa 

10 

43 

53 

0.296 

0.094 

0.140 

0.094 

kds5 

10 

86 

96 

0.288 

0.103 

0.148 

0.096 

KDS6 

10 

172 

182 

0.254 

0.093 

0.131 

0.087 

KDS7 

10 

0 

10 

0.551 

0.177 

0.313 

0.166 

kds8 

30 

0 

30 

0.240 

0.075 

0.140 

0.066 

KDSg 

0 

35 

35 

3659.1 

1343.9 

1860.3 

1195.1 

KDS10 

0 

86 

86 

2364.5 

871.9 

1189.3 

800.47 

True  Value  Parameters 

2.88 

1.09 

1.51 

1 

Table  8:  Typical  standard  errors  for  the  transition  rates  «q,  k2,  k3  and  K4 
with  a  series  of  different  types  of  data  sets. 


Figure  7,  and  also  support  our  intuitive  reasoning  about  flows  in  closed  net¬ 
works  functioning  at  capacity.  For  such  networks,  the  transition  rates  from 
one  state  to  the  other  are  obviously  highly  correlated. 

In  order  to  illustrate  the  above  theory,  we  again  used  ordinary  least 
squares  inverse  problems  with  numerous  sets  of  data.  We  performed  the 
least  square  minimization  for  several  data  sets  of  type  KDS \  with  a  total  of 
45  observations,  where  10  are  uniformly  distributed  within  the  interval  [0,  30] 
and  35  are  uniformly  distributed  within  the  interval  [40,  210]  (for  simplicity 
we  exclude  the  transition  interval  [30,40]  from  our  present  analysis).  Stan¬ 
dard  errors  for  the  estimates  of  K\,  ac2,  ft3  and  k4  in  a  typical  optimization 
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with  one  of  the  KDS\  data  sets  are  shown  in  Table  8. 

When  increasing  the  number  of  points  sampled  in  the  interval  [0,  30]  to 
30  and  keeping  the  number  of  data  points  in  [40,  210]  the  same  (the  KDS2 
type  data  set,  Table  8),  the  standard  errors  corresponding  to  each  of  the 
parameters  K\,  k2,  Kg  and  k4  decrease  considerably.  The  new  standard  errors 
range  between  61%  and  68%  of  the  standard  errors  for  the  data  set  KDS\ .  A 
decrease  in  the  standard  errors  is  also  observed  when  we  solve  the  least  square 
problem  using  data  (sets  of  type  KDS7,  KDSg  of  Table  8)  only  from  the 
interval  [0,  30],  where  the  standard  errors  obtained  with  data  set  KDSg  range 
between  40%  and  45%  of  the  standard  errors  for  KDS7.  Thus,  numerical 
calculations  fully  support  the  GSF  theory  that  increasing  the  number  of 
data  points  in  the  region  [0,  30]  results  in  more  accurate  estimates  for  the 
parameter  R. 

Next  we  increased  the  number  of  data  points  sampled  from  the  interval 
[40,210],  and  obtained  an  entirely  different  outcome.  By  comparing  the 
entries  for  the  data  sets  KDS4,  KDSg  and  KDSg  in  Table  8,  we  see  that 
there  is  little  improvement  in  the  standard  errors  when  we  successively  solve 
the  least  square  problem  with  a  fixed  number  of  10  data  points  in  the  interval 
[0,30]  and  43,  86  and  respectively  172  data  points  in  [40,210].  Also,  when 
we  attempt  to  estimate  the  parameters  using  only  data  (data  sets  KDSg  and 
KDS10 )  sampled  from  the  region  [40,210],  the  resulting  standard  errors  are 
large  (the  parameter  estimates  were  also  quite  bad  for  all  of  the  KDSg  and 
KDS\o  type  data  sets),  increasing  in  magnitude  as  the  number  of  sampled 
points  increases.  As  in  the  case  of  the  logistic  growth  model  above,  this 
is  not  surprising  if  one  recalls  the  investigations  in  [9].  It  is  an  expected 
consequence  of  introducing  redundancy  in  the  sensitivity  matrix,  which  in 
turn  makes  the  Fisher  information  matrix  ill-posed,  yielding  huge  standard 
errors. 

We  also  performed  a  similar  analysis  for  the  estimation  of  the  nodal  ca¬ 
pacities  U  and  of  the  initial  conditions  Cgk  from  data  with  Gaussian  noise 
added  as  described  above.  In  Figures  8  and  9  we  depict  the  generalized  and 
traditional  sensitivity  functions  with  respect  to  l2,  I3,  U  and  c0i,  c0 2,  c0 3  and 
C04,  respectively.  In  both  cases,  one  can  distinguish  a  small  time  interval, 
at  the  beginning  of  the  time  axis,  where  the  traditional  sensitivity  functions 
exhibit  sharp  increases/decreases  before  reaching  the  steady  state.  The  same 
time  interval  also  gives  the  region  where  the  generalized  sensitivity  functions 
with  respect  to  Cok  exhibit  sharp  increases/decreases  before  following  the 
“forced-to-one”  artifact.  Based  on  the  theory,  we  can  anticipate  that  for  the 
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Generalized  Sensitivity  Functions  with  respect  to  l2J3J4  Traditional  Sensitivity  Functions  with  respect  to  sm,  l2,lg,l4 


Figure  8:  Generalized  and  Traditional  Sensitivity  Functions  for  sm,  I2,  h,  h- 


Figure  9:  Generalized  and  Traditional  Sensitivity  Functions  for  c0 1,  c0 2,  c0 3, 

CQ4- 


estimation  of  the  initial  conditions  Co/c,  the  data  sampled  from  this  initial  time 
interval  is  the  most  informative.  Numerical  calculations  presented  in  Table  9 
confirm  this  expectation.  Indeed,  increasing  the  number  of  data  points  in 
the  interval  [0,  24]  leads  to  smaller  residual  errors  and  better  estimates  (see 
the  standard  errors  corresponding  to  CDS2  and  CDSg),  whereas  increasing 
the  number  of  data  points  in  the  interval  [34,  210]  does  not  provide  much 
improvement  (see  CDS4 ,  CDS5  and  CDSq).  Similar  to  the  findings  when 
we  estimated  the  k’s,  when  we  use  only  data  points  from  the  interval  [34,  210] 
for  the  estimation  of  c0fc  (see  CDSg  and  CDSW),  we  obtain  very  large  stan- 
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Data  points 

in 

Standard  Errors  for 

Data  Set 

[0,24] 

[34,210] 

Total 

Coi 

0)2 

0)3 

0)4 

CDS, 

8 

36 

44 

0.0075 

0.0079 

0.0077 

0.0094 

cds2 

24 

36 

60 

0.0059 

0.0064 

0.0062 

0.0088 

cds3 

8 

60 

68 

0.0083 

0.0089 

0.0079 

0.0099 

CDSa 

8 

45 

53 

0.0087 

0.0084 

0.0074 

0.0097 

cds5 

8 

89 

97 

0.0089 

0.0086 

0.0077 

0.0095 

CDS6 

8 

178 

186 

0.0082 

0.0088 

0.0076 

0.0092 

CDS7 

8 

0 

8 

0.0110 

0.0116 

0.0062 

0.0124 

CDS8 

24 

0 

24 

0.0062 

0.0071 

0.0065 

0.0090 

cds9 

0 

36 

36 

1596.50 

1434.50 

879.97 

1544.90 

CDS,  o 

0 

89 

89 

630.24 

859.35 

483.58 

1095.80 

True  Values  for  Cok 

0.0945 

0.2221 

0.6313 

0.0519 

Table  9:  Typical  standard  errors  for  the  initial  conditions  c9k  with  a  series 
of  different  types  of  data  sets. 


dard  errors  (as  well  as  unacceptable  parameter  estimates)  compared  to  the 
previous  ones. 

As  established  in  [2],  sm  is  the  parameter  to  which  the  system  (20)  is  the 
least  sensitive  overall.  In  fact  for  the  particular  values  used  for  simulations 
(sm  =  0.1039,  l2  =  0.2387,  l3  =  0.6498,  Z4  =  0.0567)  the  output  of  the 
system  (20)  does  not  depend  on  sm  at  all  (see  the  traditional  sensitivity 
functions  from  Figure  8).  From  the  particular  form  of  our  system  and  by 
practical  insight,  this  can  be  explained  intuitively  by  the  fact  that  when 
the  capacity  sm  of  the  first  node  is  too  high,  the  flow  c4  (which  replenishes 
the  network)  will  never  reach  it,  making  this  capacity  constraint  inactive 
for  our  dynamical  system.  The  consequence  is  that  for  the  given  data,  the 
entries  corresponding  to  sm  in  the  sensitivity  matrix  are  all  zero,  making  the 
Fisher  information  matrix  singular.  Intuitively,  this  simply  says  that  it  is 
impossible  to  reconstruct  a  parameter  from  data  where  the  solution  values 
do  not  depend  at  all  on  that  parameter.  This  is  the  reason  for  which  in 
Figure  8  we  present  the  generalized  sensitivity  functions  only  with  respect  to 
l2 ,  l3  and  /4.  Unlike  the  generalized  sensitivity  functions  for  n  and  Cok  which 
present  an  initial  region  with  sharp  increases/decreases  followed  by  a  region 
where  they  exhibit  a  “forced-to-one”  artifact,  the  GSF  for  l2 ,  l3  and  U  show  a 


Data  points 

in 

Standard  Errors  for 

Data  Set 

[0,24] 

[44,210] 

Total 

I2 

h 

u 

LDS1 

8 

34 

42 

0.0023 

0.0062 

0.0031 

lds2 

24 

34 

58 

0.0019 

0.0051 

0.0026 

lds3 

8 

56 

64 

0.0018 

0.0049 

0.0025 

LDS4 

8 

42 

50 

0.0020 

0.0056 

0.0028 

lds5 

8 

84 

92 

0.0014 

0.0041 

0.0020 

lds6 

8 

168 

176 

0.0011 

0.0030 

0.0015 

lds7 

8 

0 

8 

0.0071 

0.0165 

0.0086 

lds8 

24 

0 

24 

0.0032 

0.0077 

0.0039 

LDSg 

0 

34 

34 

0.0026 

0.0073 

0.0037 

LDS10 

0 

84 

84 

0.0015 

0.0044 

0.0021 

True  Values  for  l2,  l3  and  Z4 

0.2387 

0.6498 

0.0567 

Table  10:  Typical  standard  errors  for  the  capacities  l2,  I3  and  U  with  a  series 
of  different  types  of  data  sets. 


steady  increase  from  zero  to  one  throughout  their  time  courses.  We  anticipate 
therefore  that  for  the  estimation  of  the  capacities  /,  the  information  is  more 
uniformly  distributed  in  the  interval  [0,210].  Numerical  results  presented  in 
Table  10,  where  higher  number  of  data  points  result  in  lower  standard  errors, 
confirm  these  expectations. 

We  conclude  this  section  with  the  observation  that  the  analysis  and  the 
numerical  simulations  reported  here  illustrate  the  utility  of  the  GSF  in  in¬ 
vestigations  of  reasonably  complex  estimation  problems.  However,  as  in  the 
case  of  the  logistic  growth  model,  we  emphasize  that  one  should  use  care 
when  making  inferences  from  the  regions  of  monotonicity  for  these  functions. 
In  addition  to  the  regions  of  genuine  information  content,  the  GSF  often 
indicate  regions  of  steep  increase  (the  interval  [40,210]  in  Figure  7  and  the 
interval  [34,  210]  in  Figure  9)  which  are  due  simply  to  the  “forced-to-one” 
artifact  inherent  in  the  GSF  as  defined  in  [26]. 
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5  Conclusions 


From  the  analysis  presented  in  this  paper  and  in  [11,  26],  it  is  clear  that 
generalized  sensitivity  functions  are  useful  tools  in  inverse  problem  investi¬ 
gations.  The  primary  positive  feature  of  generalized  sensitivity  functions  is 
their  ability  to  suggest  regions  of  high  information  content  where,  if  addi¬ 
tional  observations  are  taken,  one  can  generally  expect  to  improve  specific 
parameter  estimates.  This  is  supported  by  our  numerical  findings  for  the 
logistic  growth  model  and  a  recently  developed  agricultural  production  net¬ 
work  model,  and  further  supports  the  initial  computational  findings  with 
examples  in  [26].  Moreover,  by  visually  investigating  the  dynamics  of  TSF 
and  GSF  curves,  we  can  potentially  identify  subsets  of  parameters  which  are 
highly  correlated.  Although  the  GSF  theory  alone  can  be  of  use  in  formula¬ 
tion  of  parameter  estimation  problems,  the  most  insight  can  be  gained  when 
the  GSF  are  used  in  conjunction  with  traditional  sensitivity  functions.  This 
will  ensure  maximum  understanding  of  parameter  sensitivity  and  can  be  of 
great  value  in  improving  data  acquisition  design. 
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